/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 *                                                                           *
 *                          Copyright (c) 2010                               *
 *                  Alen Stojanov [ ioi.alen@gmail.com ]                     *
 *                Jacobs University Bremen, Bremen, Germany                  *
 *                                                                           *
 *                          All Rights Reserved.                             *
 *                                                                           *
 *    Permission to use, copy, modify and distribute this software and       *
 *    its documentation is hereby granted, provided that both the            *
 *    copyright notice and this permission notice appear in all copies of    *
 *    the software, derivative works or modified versions, and any           *
 *    portions thereof, and that both notices appear in supporting           *
 *    documentation.                                                         *
 *                                                                           *
 *    This code is distributed in the hope that it will be useful, but       *
 *    WITHOUT ANY WARRANTY; without even the implied warranty of             *
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. THE AUTHORS       *
 *    DISCLAIM ANY LIABILITY OF ANY KIND FOR ANY DAMAGES WHATSOEVER          *
 *    RESULTING FROM THE USE OF THIS SOFTWARE.                               *
 *                                                                           *
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#include <iostream>
#include "lsp.h"

using namespace std;


LSP::LSP(QVector<TSData*> ts)
{
	this->S = ts;
}

QVector<TSPoint> LSP::calculate (int controlPointsSize)
{
	if ( projection.size() == 0 && S.size() >= 2 ) {
		/* find the control points using PAM */
		pam ( controlPointsSize );
		/* use Sammon's projection to project the points */
		sammon();
		/* use the LSP to finally calculate all projection points */
		lsp();
	}
	return projection;
}

void LSP::lsp ()
{
	/* calculate alpha values */
	TSMatrix A(S.size() + C.size(), S.size());

	/* map the position of every point to the proper control number */
	QVector<int>  S2C(S.size(), -1);
	for ( int i = 0; i < C.size(); i++) {
		int index = S.indexOf(C[i]);
		if ( index != -1 ) {
			S2C[index] = i;
		}
	}

	/* initialize the part of the matrix with control points */
	for ( int i = S.size(); i < S.size() + C.size(); i++ ) {
		/* shift the index */
		int index = i - S.size();
		for ( int j = 0; j < S.size(); j++ ) {
			if (  S2C[j] == index ) {
				/* if we have a control point, then it is one */
				A(i,j) = 1;
			} else {
				/* otherwise we have zero */
				A(i,j) = 0;
			}
		}
	}

	/* fill the L part of the matrix */
	for ( int i = 0; i < S.size(); i++ ) {
		for ( int j = 0; j < S.size(); j++) {
			if ( i == j ) {
				A(i,j) = 1;
			} else {
				A(i,j) = 0;
				if ( S2C[i] != -1 ) { /* the point is a control point */
					QVector<TSData*>& V_i = V[S[i]];
					if ( V_i.indexOf(S[j]) != -1 ) { /* p_j is in V_i */
						A(i,j) = alpha (i, j, V_i);
					}
				}
			}
		}
	}

	/* so far A has been filled, now we need to solve the system */
	TSMatrix M ( ((~A) * A).Inv() * (~A) );

	/* define two vectors one for x coordinates, one for y coordinates */
	TSMatrix bx (S.size() + C.size(), 1);
	TSMatrix by (S.size() + C.size(), 1);

	/* initialize bx and by */
	for ( int i = 0; i < S.size(); i++ ) {
		bx(i, 0) = 0;
		by(i, 0) = 0;
	}
	for ( int i = S.size(); i < S.size() + C.size(); i++ ) {
		int index = i - S.size();
		bx(i, 0) = CMDP[index].getX();
		by(i, 0) = CMDP[index].getY();
	}

	/* finally calculate the vectors */
	TSMatrix rX = M * bx;
	TSMatrix rY = M * by;

	cout << "==== Control Points ====" << endl;
	for ( int i = 0; i < CMDP.size(); i++ ) {
		cout << CMDP[i].getX() << " " << CMDP[i].getY() << endl;
	}
	cout << "========================" << endl << endl;

	cout << "======== bX bY =========" << endl;
	for ( unsigned int i = 0; i < bx.RowNo(); i++ ) {
		cout << bx(i,0) << " " << by(i,0) << endl;
	}
	cout << "========================" << endl << endl;

	cout << "========== A ===========" << endl;
	cout << A << endl;
	cout << "========================" << endl << endl;

	cout << "========== M ===========" << endl;
	cout << M << endl;
	cout << "========================" << endl << endl;

	cout << "======== rX rY =========" << endl;
	for ( unsigned int i = 0; i < rX.RowNo(); i++ ) {
		cout << rX(i,0) << " " << rY(i,0) << endl;
	}
	cout << "========================" << endl << endl;


	for ( int i = 0; i < S.size(); i++ ) {
		TSPoint tmp(rX(i, 0), rY(i,0));
		projection.push_back(tmp);
	}

}

void LSP::sammon()
{
	/* initialize the algorithm */
	for (int i = 0; i < C.size(); i++) {
		double x = qrand() / (double) RAND_MAX;
		double y = qrand() / (double) RAND_MAX;
		TSPoint tmp(x, y);
		CMDP.push_back(tmp);
	}

	TSPoint delta(0.0, 0.0);
	int maxIterations	= 100;
	double lambda		= 1;
	double d_ij			= 0.0;
	double D_ij			= 0.0;

	/* initialize the random seed */
	qsrand(0);

	for ( int iteration = 0; iteration <= maxIterations; ++iteration ) {
		cout << "Iteration: " << iteration << endl;
		/* calculate the lambda value for each iteration */
		double ratio = (double)iteration / (double)maxIterations;
		lambda = pow(0.01, ratio);
		cout << "Lambda: " << lambda << endl;
		/* create a random vectors to make sure that the values are randomly tackled */
		QVector<int> indexesI = randomPermutationVector(C.size());
		QVector<int> indexesJ = randomPermutationVector(C.size());

		/* perform the minimization */
		for (int i_itt = 0; i_itt < C.size(); i_itt++) {
			int i = indexesI[i_itt];
			for (int j_itt = 0; j_itt < C.size(); j_itt++) {
				int j = indexesJ[j_itt];
				if ( i != j ) {
					/* calculate the distance of the time series */
					d_ij = C.at(i)->distance(C.at(j));
					/* calculate the Eucledian distance of the projection points */
					D_ij = CMDP[i].distance(CMDP[j]);

					/* avoid devision with 0 */
					if (D_ij == 0) {
						D_ij = 1e-10;
					}

					/* calculate the delta correction */
					delta = (CMDP[i] - CMDP[j]) * (lambda * (d_ij - D_ij) / D_ij);

					// cout << (lambda * (d_ij - D_ij) / D_ij) << " ";

					//					cout << "i: " << CMDP[i].getX() << " " << CMDP[j].getY() << endl;
					//					cout << "j: " << CMDP[j].getX() << " " << CMDP[j].getY() << endl;
					//					cout << "---" << delta.getX() << " " << delta.getY() << endl;

					/* perform the correction */
					CMDP[i] = CMDP[i] + delta;
					CMDP[j] = CMDP[j] - delta;

					//					cout << "i: " << CMDP[i].getX() << " " << CMDP[j].getY() << endl;
					//					cout << "j: " << CMDP[j].getX() << " " << CMDP[j].getY() << endl;
					//					cout << endl;

				}
			}
		}
	}
}

void LSP::pam (int nc)
{

	/* we store the medoids in this vector */
	QVector<TSData*> medoids;
	/* we store the non-medoids in this vector */
	QVector<TSData*> non_medoids;

	/* populate the vectors, assuming the first nc points are the medoids */
	for (int i = 0; i < S.size(); i++) {
		if ( i < nc ) {
			medoids.push_back(S.at(i));
		} else {
			non_medoids.push_back(S.at(i));
		}
	}

	bool change = true;
	double cost = pamConfCost(medoids, non_medoids);

	while ( change ) {
		change = false;
		for (int i = 0; !change && i < medoids.size(); i++) {
			for (int j = 0; !change && j < non_medoids.size(); j++) {
				/* perform switch */
				TSData* tmp = medoids[i];
				medoids[i] = non_medoids[j];
				non_medoids[j] = tmp;

				/* recalculate new cost */
				double newCost = pamConfCost(medoids, non_medoids);

				/* check the new cost */
				if ( newCost < cost ) {
					/* if the cost is lower, then we are proceeding with it */
					cost = newCost;
					change = true;
				} else {
					/* otherwise we switch the medoids back */
					TSData* tmp = medoids[i];
					medoids[i] = non_medoids[j];
					non_medoids[j] = tmp;
				}
			}
		}
	}

	getNeighbourhoods(medoids, non_medoids);

	/* the medoids are actually our control points */
	C = medoids;
}

double LSP::pamConfCost(QVector<TSData*>& medoids, QVector<TSData*>& non_medoids)
{
	double tmp_distance = 0;
	double conf_cost = 0;
	for (int i = 0; i < non_medoids.size(); i++) {
		TSData* non_m = non_medoids.at(i);
		double min_distance = non_m->distance(medoids.at(0));
		for (int j = 1; j < medoids.size(); j++) {
			tmp_distance = non_m->distance(medoids.at(j));
			if (tmp_distance < min_distance) {
				min_distance = tmp_distance;
			}
		}
		conf_cost += min_distance;
	}
	return conf_cost;
}

QVector<int> LSP::randomPermutationVector(int size)
{
	QVector<int> tmp(size);
	for (int i = 0; i < size; ++i) {
		int j = floor(qrand() / (double)RAND_MAX * (size-1));
		tmp[i] = tmp[j];
		tmp[j] = i;
	}
	return tmp;
}

double LSP::alpha(int i, int j, QVector<TSData*>& V_i)
{
	double distance = 1 / ( S[i]->distance(S[j]) );
	double sum_dist = 0;

	for (int v_itt = 0; v_itt < V_i.size(); v_itt++) {
		sum_dist += 1 / (S[i]->distance(V_i[v_itt]) ) ;
	}

	return distance / sum_dist;
}

void LSP::getNeighbourhoods (QVector<TSData*>& medoids, QVector<TSData*>& non_medoids)
{
	double	tmp_distance = 0;
	V.clear();

	for (int j = 1; j < medoids.size(); j++) {
		/* create an empty vector to define each neighborhood */
		QVector<TSData*> tmp;
		/* create mapping from each control point to its neighboring set */
		V.insert(medoids.at(j), tmp);
	}

	for (int i = 0; i < non_medoids.size(); i++) {
		TSData& non_m = *non_medoids.at(i);
		double min_distance = non_m.distance(medoids.at(0));
		int min_index = 0;
		for (int j = 1; j < medoids.size(); j++) {
			tmp_distance = non_m.distance(medoids.at(j));
			min_index 	 = j;
			if (tmp_distance < min_distance) {
				min_distance = tmp_distance;
			}
		}

		V[medoids.at(min_index)].push_back(non_medoids.at(i));
	}
}

QVector<int> LSP::getConnections()
{
	QVector<int> connections(S.size(), -1);
	if (V.size() != 0) {
		QMap <TSData*, QVector<TSData*> >::iterator it;
		for (it = V.begin(); it != V.end(); ++it) {
			int a = S.indexOf(it.key());
			connections[a] = a;
			for (int i = 0; i < it.value().size(); ++i) {
				int b = S.indexOf(it.value().at(i));
				connections[b] = a;
			}
		}
	}
	return connections;
}

LSP::~LSP()
{
	// Do nothing
}
